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SUMMARY 


We consider the numerical solution of the single-group, steady state, isotropic transport 
equation. An analysis by means of the moment equations shows that a discrete ordinates Sn 
discretization in direction (angle) with a least squares finite element discretization in space does not 
behave properly in the diffusion limit. A scaling of the S N equations is introduced so that the least 
squares discretization has the correct diffusion limit. For the resulting discrete system a full 
multigrid algorithm was developed. 


1. INTRODUCTION 


The single-group, steady state, isotropic form of the Boltzmann transport equation for one 
dimensional slab geometry is given by ( Lewis and Miller [6] ) 


+ #*) “ \ / x ’ ^ 

< -1 
= gM for p> 0 

k Tp(b,n) = 92 (v) fo r ^ < 0 


q{x,n) I 


(l.i) 


where x € [a, b } and p 6 [-1, 1]. When * t -» oo and % - 1, which is the so called diffusion hml, 
this equation becomes singular. The limit operator (/ - P), where P denotes the operator 
Pip(x,n) = \ f\ ip{x, p')dp', has in its nullspace all functions that are independent of angle p. 

Moreover' in this limit transport theory transitions into diffusion theory in the following way. 

Let e be a small parameter. Substituting a t by L, a a by (^ — £<r a ), where o a is 0(1), and scaling the 
right hand side by e, equation (1.1) becomes 

pA + -I - (--£<Ta) p|t/>(x,p) = eq{x). (1-2) 

In addition, it is assumed that the external source q is independent of p. As a consequence of this 
parameterization the diffusion limit is now equivalent to the limit £ -» 0. By expanding the solution 
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of (1.2) as ^ ^ . 

ip{x,n) = +f^ e 3 tp j (x,^) 

i=i 


it can be shown ( Larsen [2] [3], Larsen, Morel and Miller [4] ) that some mean-free-paths away 
from the boundary the zeroth order term, i/j°(x, /z), is independent of /z and is a solution of the 
following diffusion equation 


1 d?<fi(x) 
3 ch? 


+ cr a <f>{x) = q(x). 


(1.3) 


Thus, in the diffusion limit the solution of the transport equation will converge to the solution of a 
diffusion equation. 

For the numerical solution of (1.1) it is important to find a discretization that has the same 
property; i.e. , for diffusive regimes ( a t large, ^ ~ 1 ) the difference scheme for the transport 
equation must approximate a diffusion operator. 

In the last two decades a large amount of work was dedicated to developing special 
discretizations for the transport equation that have the correct behavior in the diffusion limit. 
Among them are the Diamond Difference scheme [6], the Linear Discontinuous scheme [1], and the 
Modified Linear Discontinuous scheme [5]. In the one-dimensional case their implementation is 
straightforward, but their extension to higher dimensions is difficult. 

In this paper we try to develop a general framework for finding discretizations for the transport 
equation that have the correct behavior in the diffusion limit. In Section 2 we describe the discrete 
ordinates Sn discretization in angle and a least squares finite element discretization in space and 
discuss why this simple approach does not behave properly in the diffusion limit. A scaling 
technique for the transport equation is introduced in Section 3 that yields a least squares 
discretization with the proper diffusion limit. In Section 4 we present numerical results based on a 
full multigrid solver for the resulting discrete system. In Section 5 we draw conclusions and suggest 
further applications of the scaling technique. 


2. DISCRETIZATION 

For the discretization in angle we use the standard discrete ordinates Sn method. In the case of 
one-dimensional Slab geometry, this is a Galerkin discretization with normalized Legendre 
polynomials as basis. That means we are looking for a flux solution that has an expansion in the 
first N normalized Legendre polynomials, 

N - 1 

( 2 . 1 ) 

1=0 

Since the normalized Legendre polynomials form an orthonormal basis for 

■^ 2 ([— 1, 1]), the moment coefficients 4>i are given by the following integral, which can be 




A *-* 


,T'i *< 
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written exactly as a sum by using a Gauss quadrature rule. We have 

M x ) = \J 

( 2 . 2 ) 

N 

i - 1 

Here Hj denotes the Gauss quadrature points and denotes the Gauss quadrature weights. 

By introducing the vector notation 

^ = (ip(x, Hi ), . . . , ip(x, Mw)) T > = (0o( x )? • ■ ■ i^N-ii 30 )) 


and defining matrices T and Q. as 

\T\ij = Pi- lOij); n = diag(wi, ..., 0 J N ) 

the relationship (2.1) and (2.2) between the flux £ and the moments $ can be written as 

$ - TQ.% ( 2 - 3 ) 

y _ j-t$ (2.4) 


As a result of the Galerkin discretization of (1.1) with the ansatz (2.1) we obtain the S N 
equations 


( Mi 


Lty = e 




— + 7* - (1 - eV„) Hi = £ 2 £, 

OX 


(2.5) 


Hn 


where 


fl = (l,...,l) T (wi f ... f WAr). 

When we insert (2.4) into (2.5) and multiply by TO from the left, we get the moment equations 



e 2 (J a 

£b °T* 

0 

0 


M$ = 

£b o £ 
0 

1 

£b '-§i 

£bl 'dx 

1 

• o 



. 




- 


= e 2 q 


(2.6) 


with 


b i = 


3 + 1 


+ !) 2 - 1 

Normally, the computations are done in the flux representation (2.5) since in this representation 
the boundary conditions are equal to simple Dirichlet boundary conditions. However, as we will see 
later, the moment equations are very useful for theoretical insight. In the following the flux 
operator is denoted by L and the moment operator by M. 
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For the spatial discretization of the Sn equations we use a least squares finite-element method 
based on the functional 

b 

F ($L) = J - ?) , dx, (2.7) 

a 

and piecewise linear continuous elements rjk as basis functions for each component of 2.. In the 
following S h denotes the finite dimensional space S h = span {r/ fc } and (S h ) N denotes the space of 
iV-tuples whose elements are in S h . 

The advantage of this approach is that a least squares discretization converts the Sn equations, 
which are a coupled system of first order equations, into a self-adjoint variational formulation 
Based on this variational formulation Multi-Level Projection Methods [7] can be applied in order to 
guide the development of a multigrid solver for the resulting discrete system. 

Unfortunately, this discretization does not behave correctly in the diffusion limit. In order to 
explain this fact, we use the moment equations, since fa = ip in the diffusion limit. Further, it is 
easy to see by using the relationship (2.3), (2.4) and the identity T T Tfl = I that 

6 

min / (ft ( L ^ - g) , D J' — q) ., dx <=$■ 

*€(S h )X J ' V — -V — 2-/ m* 


mm 

*s(S h ) 


o 

in 


which justifies why it is also possible to look at the least squares discretization of the moment 
equations. 

In the S 2 case with cr a = 0, for example, a least squares discretization of the moment equations 
results in the following discrete system 





+ (v,v) 


(*) = (>,*>)’ 


where, for example, ( 77 , rj) is a mass matrix and ( 77 ', rj) a stiffnes matrix with elements 

_ ) dr)k(x) 


[M] m = / *?*<*)!»(*)&:, [(v',v)] k ,i = / ^^m(x)dx 


and 


(V.<7o)= ( ■ • ■ , J dV ^ X \ Q (x)dx, . 


dx 

Forming the Schur complement we get the following equation for <p$ 

-2 




(v,v) + j(v',v') 


(*7>V)J 00 




(WiV) + 


(v'iQo)- 


( 2 . 8 ) 
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For sufficient small e we have 

+ ={v,v)~ 1 + o( £A )- ( 2 - 9 ) 

Plugging (2.9) into (2.8) and dividing by e 2 leads to 

o ' — — v — 

(*i) 

+^-(*?', 7 ?) (t?, i ?) -1 (*?\ »?')(»?» f?)" 1 V) + 0(e 4 ) } <Po ( 2 - 10 ) 

= -y(*/.*7)(*?.»7)" 1 (»? , i«9) + 0 ( e4 )' 

In the limit e -» 0, the solution approaches a valid discretization for the diffusion equation (1.3) 
only if the term (*1) vanishes identically. For piecewise linear, continuous basis elements this term 
does not cancel out and becomes the leading term in the equation. Consequently, in the diffusion 
limit (2.10) is an approximation for <%' = 0, which results in a linear solution, connecting the 
boundary conditions. In general, the term (*1) does not have the proper behavior unless the mass 
matrix, (17,17), is lumped, that is, replaced by a diagonal matrix. 

3. SCALING 

A closer look at the moment equations (2.6) shows that this system is unbalanced. There are 
0(e 2 ), 0(e) entries as well as 0(1) entries. The idea is to scale this system before the discretization. 

First, let us consider the case a a ± 0. In our inner product the adjoint moment operator, when 
homogeneous boundary conditions are assumed, is given by 

M* = 


Scaling the moment equations by 

S = 


and forming the normal equations results in 

M*SM§> = e 2 M*Sq <==> 
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/ F 3 rr ^ 

0 

-^bobi-gjr 

0 

... \ 

0 

E <r a dl* £ °1 dx* 

0 

~^bib 2 ^s: 


-e 3 bobi 

0 

e-e 3 (6f +b 2 2 )£ I 

0 


0 

—e 3 b 

0 

e - £ 3 (&2 + 




( £% 


• / 


z gl a i\ 
vJ cr a 


[ ’ J 


Note that (3.2) already has the correct limit equations for the moments on its diagonal and all first 
order derivatives are eliminated. 


Applying a Galerkin discretization to (3.2) and forming the Schur complement leads after 
division by e 3 to the following discrete equation for 


{<7a(r?, rj) + ifo', rf) + 0(e 2 )} = (r/, q 0 ) + 0(e 2 ), 

which is a valid discretization of the corresponding diffusion equation (1.3) in the limit e — > 0. 

When we define b k j = €jVk(x), where e, denotes the j-th canonical unit vector of IR N , we can write 
the Galerkin discretization of (3.2) as follows 


b 

J {M'S (MJ>- qj dx = 0. (3.3) 

a 

Assuming homogeneous boundary conditions and splitting S = y/SVS, (3.3) is equivalent to 

6 

V&M J (ys (Mi - j) , VSMb kJ ) nK dx = 0 


•*=*• (Mi - j) , x/S (Mi - $)) dx, 

a 

which is a least squares discretization of the moment equations, scaled by VS. Consequently, a 

least squares discretization of the moment equations, scaled by </S, also has the correct behavior in 
the diffusion limit. 

Notice that (3.2) has the proper behavior for any cr a ^ 0. The second equation contains J- on 
both sides and yields the proper solution for ^ as a a — + 0. 


398 


On the other hand, if cr a = 0 the scaling (3.1) cannot be applied. In this case, scaling the 
moment equations by 

/ 1 \ 


S 0 = 


a 


a 


\ 


with 


a = e p , where p > 2 


(3.4) 


(3.5) 


and forming the normal equations results in 


-EOLb 2 ^ 
e :aO 0 dP 

-Eab 0 £ 

— e 2 a:6o&i 

0 


eb o£ 

a-E 2 {bl + ab\)¥p_ 

0 

-e 2 ab 162^5 


-e 2 abobi-§p 

0 a - 

-A »(* 5 + 

0 


0 

—E 2 ab 1&2 

0 

a - £ 2 a{bl + bj) £? 

'■ J 


( 0 \ 

-^*0^ 

0 


V J 

A Galerkin discretization of (3.6) leads to the following discrete system 


(3.6) 


where 


and An = 


( 0 \ 


E 2 abl{r)',r]') 

Am W - 

£ 3 b 0 (r)’, q 0 ) 

Aio 

An)- ~ 

0 


V J 


A 0 i = (- eabo{r),r]’),e 2 ab 0 bi(r) , ,r)'),0 , . . .) 
Aio = (Eabo(r),r) , ),E 2 ab 0 bi(rj',T} , ),0,.. ,) T 


(3.7) 


^ oc{r} i rj)+e 2 {bl+ab t 2 )(r]' 1 rj f ) 0 e 2 ah b 2 (rf / , ij') 0 . ^ 

0 <2(77, v ) + f 2 a(6| + rf ) 0 { w r , v ') 

e 2 ab\ b>(T}‘, 77') 0 a(f?, r/) + £ 2 q( 6J + 6?) 0 

0 c 2 abib>{q , 9 ri r ) 0 a(r} t rj) + e 2 a(i>? 4 - 6§) 
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■J 


For sufficient small e we have 
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with Q\ — 0(~r) and Q 2 = 0(^-). Using this expansion when forming the Schur complement in 
(3.7) we get 

( etboirf.qo) \ 

\e 2 abl{r}' ,T)') - <f>o = -^oi^n 




) 


which after some algebra becomes 


{e 2 abl{ri',ri') - a 2 {T],p){r]' ,T)') *( 77 , t/) + 0(e 4 a)} $ 
= -£ 2 a(7 7 ,r 7 ')(r ? ',r 7 ')- 1 (V,go)+0(a 2 ) 


(3.8) 


Because of (3.5) we have ^ — » 0, so that (3.8) is a valid discretization of the corresponding 
diffusion equation in the diffusion limit. Using the same argument as above it follows that the least 
squares discretization of the moment equations, scaled by y/So, also has the correct behavior in the 
diffusion limi t _ „ 

As mentioned before, the computations are done in the flux representation. Therefore, we need 
to transfer the scaling of the moment equations to a scaling for the Sn equations. By means of the 
relationship (2.3) and (2.4) it follows that a least squares discretization of the moment equations, 
scaled by \/S is equivalent to a least squares discretization of the Sn equations, scaled by 

T t VSTQ = — L -.R + y/i(I-R). 


A further multiplication by yfeu~ a leads to the following scaling in the case cr a 7 ^ 0 

R + £y/G~a (I ~ R) • 


(3.9) 


Similarly, in the case cr a = 0 the least squares discretization of the moment equations, scaled by 
y/So with p = 4, is equivalent to a least squares discretization of the Sn equations, scaled by 


R + e 2 ( I-R ). 


(3.10) 


In order to avoid an “if else” in the computations it is possible to combine the scalings (3.9) and 
(3.10) to 

R+ (e^ + e 2 ) (I -it). (3.11) 


4. NUMERICAL RESULTS 


For the solution of the discrete system that results from a least squares discretization of the Sn 
equations scaled by (3.11), a full multigrid in space algorithm with 

• standard coarsening in space by doubling the mesh width, 

• /i-line red-black smoothing, 
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Scalar Flux: ( mesh size h = 1 .25 ) 



Figure 4.1: Scalar flux solution of problem (4.1) 

• full weighting, 

• linear interpolation, 

was developed. The full multigrid process starts by solving the problem on the coarsest level and 
uses this solution as a starting guess for the next finer level, where a single V-cycle is performed. 
Recursively, the solution process proceeds from coarser levels to finer levels by halving each grid 
cell, using the coarse level solution as a starting guess and performing a single V-cycle on the next 
finer mesh. This algorithm yields V-cycle convergence rates that are below 0.09. Therefore, by 
performing one full multigrid V-cycle, a solution with an error on the order of the truncation error 
is obtained (cf. [7]). 

As test problem we used the same problem that was used by Larsen, Morel and Miller in [4], 
which is shown below: 

/i ,~ ' l + 100 -07 - 100 = 0.01 

< dx «-i l. (4.1) 

ipj (0) = 0 for > 0 

k ^j(IO) = 0 for fij < 0 , 

In our parametrization (1.2) this implies e = 0.01, <r a = 0 and q = 1.0. The exact solution of the 
corresponding diffusion equation is 

4>{x) = - -x 2 + 15x, 
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Scalar Flux: ( mesh size h = 1.25 ) 



which is plotted in Figure 4.1 in solid. The least squares solution of the scaled S N equations, 
computed by one full multigrid V-cycle, is shown in Figure 4.1 by the crosses. We see that this is a 
very satisfactory result, especially when we take into consideration that we used a mesh size of 1.25, 
which is much larger than e. 

Finally, we mention that the least squares discretization of the S N equations without scaling will 
give the zero solution for problem (4.1), indicated by the stars in Figure 4.1. 

For the sake of completeness we present in Figure 4.2 the results for the test problem 

I A*J^ + 100 iij - 99.99 ]T = 0.01 (l + 1 5x) ) 

■ «•*> 

ipj( 10) = 0 for Hj < 0 J 

where a a = 1.0, e — 0.01, q = 1 — |x 2 + 15a:. The exact solution of the corresponding diffusion 
equation is the same as for problem (4.1) and is again plotted in Figure 4.2 in solid. The least 
squares solution, computed by 1 full multigrid V-Cycle is given in Figure 4.2 by the crosses and the 
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solution for the least squares discretization of the Sn equations without scaling is given by the stars. 

5. CONCLUSIONS 

From the analysis in Section 3 and the numerical results presented in Section 4, we conclude that 
the least squares discretization of the scaled Sn equations has the proper behavior in the diffusion 
limit. 

Further, we point out that the scaling can be used in the case of nonhomogeneous material, 
where < 7 t or, equivalently, e are discontinuous, because the equations are only scaled from the left 
side so that no derivatives are applied to the scaling operator. 

Adaptive refinement can be combined with the full multigrid solver in a natural way. Areas of 
new refinement can be identified by examining the difference in the solution for two consecutive 
grids. This is especially important for nonhomogeneous material, where interior layers may exist. 

Numerical results show that with a slightly different scaling both a Galerkin finite element 
formulation with piecewise linear elements and an Upwind Difference discretization of the Sn 
equations also have the correct diffusion limit. We believe that this scaling approach will result in a 
general framework for the development of discretizations that posses the correct diffusion limit. 

Finally, we hope to apply the scaling techniques developed here to higher dimensional problems. 
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